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Abstract 

Deposition  of  carbon  on  conventional  anode  catalysts  and  formation  of  large  temperature  gradients  along  the  cell  are  the  main  barriers  for 
implementing  internal  reforming  in  solid  oxide  fuel  cell  (SOFC)  systems.  Mathematical  modeling  is  an  essential  tool  to  evaluate  the  effectiveness 
of  the  strategies  to  overcome  these  problems.  In  the  present  work,  a  three-dimensional  model  for  a  planar  internal  reforming  SOFC  is  developed. 
A  co-flow  system  with  no  pre-reforming,  methane  fuel  utilization  of  75%,  voltage  of  0.7  V  and  current  density  of  0.65  A  cm-2  was  used  as  the 
base  case.  The  distributions  of  both  temperature  and  gas  composition  through  the  gas  channels  and  PEN  (positive  electrode/electrolyte/negative 
electrode)  structure  were  studied  using  the  developed  model.  The  results  identified  the  most  susceptible  areas  for  carbon  formation  and  thermal 
stress  according  to  the  methane  to  steam  ratio  and  temperature  gradients,  respectively.  The  effects  of  changing  the  inlet  gas  composition  through 
recycling  were  also  investigated.  Recycling  of  the  anode  exhaust  gas,  at  an  optimum  level  of  60%  for  the  conditions  studied,  has  the  potential  to 
significantly  decrease  the  temperature  gradients  and  reduce  the  carbon  formation  at  the  anode,  while  maintaining  a  high  current  density. 

©  2007  Elsevier  B.V.  All  rights  reserved. 

Keywords:  SOFC;  Direct  internal  reforming;  Methane;  Anode-supported;  Recycling 


1.  Introduction 

Solid  oxide  fuel  cells  (SOFC)  have  attracted  considerable 
interest  during  the  past  decade  as  highly  effective  sources  of 
electrical  energy  [1],  This  type  of  fuel  cell  operates  at  rela¬ 
tively  high  temperatures  (650-950  °C)  with  hydrocarbons  or 
hydrogen  as  the  fuel  at  the  anode  and  air  as  the  oxidant  at  the 
cathode.  Internal  reforming  of  methane  within  an  SOFC  to  pro¬ 
duce  hydrogen  in  situ  can  maximize  the  efficiency  and  simplicity 
of  the  process  by  combining  the  reforming  and  electrochemical 
processes.  The  two  main  drawbacks  to  internal  reforming  are  the 
large  temperature  gradients  generated  within  the  fuel  cell  and  the 
deposition  of  carbon  on  the  anode  that  leads  to  deactivation  of 
the  electrocatalysts. 

Steam  reforming  is  a  highly  endothermic  reaction,  whereas 
electrochemical  oxidation  is  an  exothermic  reaction  [2-5],  At 
the  inlet  of  a  cell  where  the  fuel,  typically  methane,  concen¬ 
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tration  is  high,  the  reforming  rate  dominates  such  that  the  inlet 
region  cools  below  the  feed  temperature.  As  the  fuel  flows  into 
the  cell  the  relative  amount  of  hydrogen  to  methane  increases 
and  the  electrochemical  reaction  eventually  dominates  result¬ 
ing  in  a  temperature  increase  downstream  in  the  cell.  The  steep 
temperature  gradients  that  result  from  the  localized  heating  and 
cooling  within  the  cell  can  lead  to  significant  thermal  stresses  in 
the  solid  structure  and  potentially  to  a  system  failure  from  crack 
formation  in  the  cell. 

The  second  drawback  of  steam  reforming  on  conventional 
nickel/yttria-stabilized  zirconia  (Ni/YSZ)  anodes  is  the  tendency 
of  Ni  to  catalyze  carbon  formation  [6].  Eventual  deactivation  of 
the  catalyst  occurs  if  the  rate  of  carbon  deposition,  occurring 
through  reaction  (1),  exceeds  the  rate  of  its  removal,  occurring 
through  reactions  (2)  and  (3)  [7-10]: 


CH4  -»  C  +  2H2 

(1) 

c  +  h2o  ->  CO  +  h2 

(2) 

C  +  O2-  -*  CO  +  2e“ 

(3) 
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Nomenclature 

A 

area 

CP 

heat  capacity  at  constant  pressure 

D 

diffusion  coefficient 

Dh 

hydraulic  diameter 

E 

energy 

£° 

reference  reaction  potential 

Eocw 

open  circuit  potential 

F 

Faraday  no. 

G 

heat  transfer  parameter 

h 

convective  heat  transfer  coefficient 

H 

height 

i 

current  density 

k 

exchange  current  density 

k 

conductive  heat  transfer  coefficient 

L 

cell  length 

«e 

no.  of  electrons  transferred  in  the  reaction 

N 

no.  of  moles 

Ni 

molar  flux  of  species  i 

Nu 

Nusselt  no. 

P 

pressure 

Prandtl 

Prandtl  no. 

R 

universal  gas  constant 

R5,6J 

reaction  rate  (molm-2  s-1) 

Re 

Reynolds  no. 

t 

thickness 

T 

temperature 

u 

linear  velocity 

w 

width 

X 

mole  fraction 

Z 

distance  in  ^-direction 

Greek  letters 

Y 

heat  transfer  parameter 

0 

potential  (voltage)  drop 

vij 

stoichiometric  coefficient 

Subscripts 

a 

anode 

act 

activation 

air 

air 

c 

cathode 

cb 

catalyst  bed 

cone 

concentration 

e 

electrolyte 

f 

fuel 

fc 

fuel  channel 

ic 

interconnect 

s 

solid 

1 

total 

Superscripts 

Eff 

effective 

Ref 

reference 

TPB 

triple  phase  boundary 

At  lower  temperatures  -  typically  below  973  K  -  dispro¬ 
portionation  of  carbon  monoxide  (the  so-called  Boudouard 
reaction)  can  also  serve  as  a  carbon  forming  reaction  as  will 
be  discussed  later. 

The  rate  of  carbon  formation  is  affected  by  the  steam  to  hydro¬ 
carbon  or  carbon  ratio  [8,11,12],  the  current  density  [9,13,14], 
and  the  hydrogen  to  hydrocarbon  ratio  [15].  Higher  steam  to 
carbon  ratios  and  higher  current  densities  decrease  the  rate  of 
carbon  formation  by  enhancing  removal  of  carbon  from  the 
surface  through  reactions  (2)  and  (3)  [9,13],  Higher  current 
densities  result  in  more  oxygen  anions  diffusing  through  the 
electrolyte  to  the  anode  structure  and  production  of  more  steam 
within  the  pores  of  the  anode  through  the  electrochemical  reac¬ 
tion.  Diffusion  of  oxygen  anions  is  limited  to  the  active  region 
of  the  anode  [16]  and  will  not  reduce  carbon  formation  at  the 
electrochemically  inactive  regions  of  thicker  anodes  in  anode- 
supported  systems.  A  higher  ratio  of  hydrogen  to  methane  will 
also  decrease  the  rate  of  carbon  formation  by  promoting  the 
reverse  of  reaction  (1)  [17,18], 

A  promising  way  to  resolve  the  problems  of  both  carbon  for¬ 
mation  and  steep  temperature  gradients  in  SOFC  is  to  recycle 
part  of  the  outlet  stream  from  the  anode  [19-21],  This  exhaust 
stream  contains  hydrogen  (H2,  ~15%),  water  and  a  relatively 
high  carbon  dioxide  (CO2)  to  carbon  monoxide  (CO)  ratio. 
Recycling  of  the  outlet  anode  stream  adjusts  the  inlet  feed  com¬ 
position.  Specifically,  the  concentration  of  methane  (CH4)  is 
decreased  and  the  concentrations  of  H2  and  CO2  are  increased. 
These  changes  result  in  a  higher  current  density  at  the  inlet, 
lower  rate  of  carbon  formation  [2],  and  decreased  steam  make¬ 
up  requirement  [19].  In  addition,  temperature  gradients  in  the 
cell  are  reduced  because  of  the  increase  in  the  relative  rate  of 
the  electrochemical  reaction  (reaction  (3))  versus  the  rate  of  the 
reforming  reaction  (reaction  (5),  below).  The  amount  of  car¬ 
bon  formation  is  reduced  from  the  higher  concentrations  of  both 
hydrogen  and  carbon  dioxide  in  the  feed.  Hydrogen  impedes 
reaction  (1),  while  carbon  dioxide  (CO2)  promotes  the  reverse 
Boudouard  reaction  [22,23]  as  shown  in  reaction  (4): 

C02  +  C  — >  2CO  (4) 

In  order  to  better  understand  the  operation  of  SOFC,  a  num¬ 
ber  of  theoretical  models  have  been  developed.  Some  studies 
have  modeled  various  SOFC  geometries  for  electrode-  [1,24] 
and  anode- supported  [25-27]  cells  with  either  internal  [25,26]  or 
external  reforming  [27],  The  results  of  these  studies  indicate  the 
formation  of  large  temperature  gradients  along  the  length  of  the 
cell  when  internal  reforming  is  employed.  The  models  published 
thus  far,  however,  contain  deficiencies  such  as  modeling  in  only 
one  or  two  dimensions  rather  than  in  three  dimensions,  assum¬ 
ing  reactions  occurring  in  the  gas-phase  are  representative  of  the 
solid-phase  reactions,  and  neglecting  the  effect  of  diffusion  of 
gases  through  the  PEN  (positive  electrode/electrolyte/negative 
electrode)  structure  in  all  directions.  These  deficiencies  affect 
the  reliability  of  the  results  and  limit  their  usefulness. 

In  this  study,  we  have  developed  a  rigorous  model  of  an  SOFC 
with  direct  internal  reforming  in  order  to  investigate  the  prob¬ 
lems  of  carbon  formation  and  temperature  gradients  within  a  cell 
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Fig.  1.  Schematic  of  modeled  system:  (a)  3D  view  of  three  cells  in  a  stack  and 
(b)  front  view  of  half  of  one  unit  cell.  The  letters  in  (b)  are  those  referred  to  in 
the  boundary  conditions  in  Table  2. 

and  stack.  The  three-dimensional  model  includes  electrochem¬ 
ical  reactions,  and  heat  and  mass  transfer  between  the  solid  and 
gas  phases,  and  is  solved  using  COMSOL  Multiphysics  3.2.  The 
composition  of  the  feed  gas  mixture,  in  particular  the  concentra¬ 
tions  of  H2,  CH4,  H2O  and  CO2,  is  the  main  parameter  varied  in 
this  study.  Our  results  demonstrate  that  recycling  of  the  anode 
exhaust  gases  is  indeed  a  promising  technique  for  regulating  the 
temperature  gradients  and  carbon  formation  within  the  cell. 

2.  System  and  model  development 

A  three-dimensional  mass  and  heat  transfer  model  was  devel¬ 
oped  for  an  anode- supported  cell  of  planar  geometry,  as  shown 
in  Fig.  1.  The  PEN  structure  is  between  interconnects  that  form 
channels  for  air  and  fuel  flow.  In  the  model,  the  air  and  fuel 
flow  co-currently  over  opposite  sides  of  the  PEN  structure.  The 


Table  1 

System  specifications 

Gas  inlet  temperature 
Current  density  (i) 

Air  to  fuel  ratio 
Anode  thickness  (4) 

Cathode  thickness  (4) 

Flow  arrangement 
Operating  voltage 
Cell  dimensions 
Fuel  utilization 
Electrolyte  thickness  (4) 
Channel  dimensions  (tufc  x  Hfc) 
Steam  to  carbon  ratio  in  fuel 


1073  K 
0.65  A  cm"2 

1000  p.m 
50  p.m 

Co-current  flow 
0.7  V 

10  cm  x  10  cm 
75% 

10p,m 
5  cm  x  2  cm 


choice  of  co-current  flow  is  discussed  in  Section  3.  Calculations 
were  performed  on  one  unit  cell  in  the  stack.  It  will  be  shown 
that  one  unit  cell  is  representative  of  the  behavior  of  the  entire 
stack,  except  at  the  edges  of  the  stack  where  the  boundary  con¬ 
ditions  may  be  different.  As  shown  in  Fig.  la,  there  is  a  plane 
of  symmetry  in  the  centre  of  each  unit  cell  and  so  the  calcula¬ 
tions  required  can  be  further  reduced.  In  the  presented  results, 
the  x,  y  and  z  dimensions  have  been  normalized  between  0  and 
1,  unless  otherwise  noted.  A  summary  of  system  specifications 
is  presented  in  Table  1. 

2.1.  Mass  balances 


The  following  three  surface  reactions  were  considered  in  the 
model: 


CH4  +  H2OrT-5'eCO  +  3H2 

(5) 

_  reversible _ 

co  +  h2o  co2  +  h2 

(6) 

H2  +  02-^  H20  +  2e" 

(7) 

The  steam  reforming  and  water-gas  shift  reactions  (5)  and 
(6),  respectively,  generate  H2  while  the  electrochemical  oxida¬ 
tion  reaction  (7),  consumes  H2  and  produces  steam  at  the  anode. 
As  the  rate  of  CO  oxidation  is  slower  than  the  rate  of  H2  oxi¬ 
dation  [27-29]  and  the  concentration  of  CO  is  lower  than  the 
concentration  of  H2,  electrochemical  oxidation  of  CO  was  not 
included  in  the  model. 

Table  2  summarizes  the  equations  and  boundary  conditions 
for  the  mass  transfer  calculations.  Eq.  (8)  describes  diffusion 
within  the  anode  structure.  All  reactions  occurring  in  this  system 
are  assumed  to  take  place  on  the  solid  surface  (i.e.,  not  in  the  gas 
phase).  The  terms  on  the  left-hand  side  of  the  equation  account 
for  the  diffusion  in  the  solid  phase  while  the  terms  on  the  right- 
hand  side  account  for  the  change  in  the  number  of  moles  of  each 
species  through  the  relevant  reactions.  The  last  term  accounts 
for  the  change  in  the  total  number  of  moles  from  the  steam 
reforming  reaction,  Eq.  (5),  in  which  the  ratio  of  the  moles  of 
products  to  the  moles  of  reactants  is  two.  The  solid  structure  is 
assumed  to  be  insulated  from  all  sides  and  mass  transfer  is  only 
permitted  from  the  anode  to  the  anode  gas  channel  as  shown  in 
the  boundary  equations  in  Eq.  (9). 

Eq.  (10)  describes  convective  mass  transfer  in  the  anode  gas 
phase.  Diffusion  of  the  gases  through  and  from  the  solid  structure 
is  accounted  for  in  this  equation  and  the  boundary  conditions 
given  in  Eq.  (11).  Diffusion  in  the  direction  of  flow  is  neglected 
and  the  only  boundary  condition  in  that  direction  is  the  initial 
concentration  of  each  species. 

Multi-component  diffusion  coefficients  were  used  for  mass 
transfer  in  the  gas  phase,  while  Knudsen  diffusion  was  used  for 
the  diffusion  through  the  porous  solid  structure  of  the  anode 
[30],  The  effective  diffusion  coefficient  (D?ff)  is  a  function  of 
the  tortuosity  of  the  solid  and  the  gas-solid  interactions  as  given 
in  Eq.  (13): 

1  _  i  1 

£jeff  —  £jKnudsen  + 


A',gas£/l 


(13) 
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Table  2 

Equations  and  boundary  conditions  used  for  mass  transfer  calculations 


Within  the  anode  structure 

Df(d2Pi/dx2  +  SPPi/dy2  +  d2Pt/dz2) 


-  £ 

.7=5,6, 7 
accordingly 


or;  or;  o  r;  or;  or; 

—  =0  —  =  0  — -  =  —Dio — -  —  =0 


Gas  phase  (anode  channel) 
1  /  d(nt(Pi/Pt))\ 


s  dy  3 z 

(d2Pi/dx2  +  a2 Pi  +  ay2)  _ 


PCH4=0.33Pt  Ph2o  =  0.67P,  Aj2  =  P 


Pi(2Ri) 


(8) 


(9) 

(10) 

(11) 

(12) 


where  /j!<nudsen  is  the  Knudsen  diffusion  coefficient  for  each 
species  i.  A, gas  the  gas  phase  diffusion  coefficient  of  species  i, 
e  the  porosity,  and  x  is  the  tortuosity. 

The  steam  reforming  reaction  rate  was  determined  using  a 
simple  reaction  model  of  the  form  shown  in  Eq.  (14)  [31]: 

Rsr  =  &o*ch4  Pt  eTEm/RTa'  (14) 

In  this  equation,  the  rate  of  reaction,  I?sr,  depends  on  the 
frequency  factor,  ko,  the  partial  pressure  of  methane,  xqw4-Pu 
activation  energy  of  the  reaction,  Act,  and  temperature  of  the 
catalyst  bed,  Tc b.  The  water-gas-shift  (WGS)  reaction  is  typi¬ 
cally  assumed  to  be  at  equilibrium  throughout  the  cell  [26,32], 
In  the  present  work,  however,  simple  kinetic  models  were  used 
for  both  forward  and  backward  reactions,  as  suggested  by 
Lehnert  et  al.  [26], 


2.2.  Heat  balances 

The  unit  cell  is  assumed  to  be  thermally  isolated  and  thus 
effectively  insulated  from  all  sides.  Heat  is  generated  or  with¬ 
drawn  from  reactions  occurring  on  the  surface  and  within  the 
solid  catalyst  structure.  The  reactions  considered  do  not  occur 
in  the  gas  phase  at  the  temperature  of  interest  for  this  study 
(i.e.,  1073  K).  Heat  is  exchanged  between  the  solid  and  the  gas 
streams  by  convection  and  through  transport  of  reactants  and 
products  from  or  through  the  solids.  Across  the  solids,  heat  con¬ 
duction  was  considered  dominant.  The  equations  and  boundary 
conditions  for  the  heat  transfer  calculations  are  given  in  Table  3. 

The  electrochemical  oxidation  consists  of  two  half-cell  reac¬ 
tions  occurring  separately  in  the  anode  and  cathode,  and  heat 
is  generated  within  the  electrolyte  and  the  electrodes.  The  heat 
generated,  however,  is  assumed  to  be  produced  within  the  anode 
structure  because  the  cathode  and  electrolyte  are  very  thin,  rel¬ 


iable  3 

Equations  and  boundary  conditions  used  for  heat  transfer  calculations 


Within  the  PEN  st 


fa  2ts  a2Ts  a2rs) 


x  =  0,  Wfc  z  =  0,  L 


k*Y=  ArAnr  -  Ts)  +  JVq jCy.o^air 
Within  the  interconnects  structure 


7  S2  Tic 


=  0 


a  ts  '  -» 

-ks  =  hf(Tf  -  Ts)  -  NiCpjTi 


'1 f 

x  =  0,  Wfc  z  =  0,L  At :  ab,  ef 

ar  „  3T  a t,c  ,  ,  aTic  ,  ars 

-  =  0  -  =  0  kic—  =  /!air(7iir  -  7/c)  kic  —  =  -ks  — 

dx  a z  3x  3y  dy 

Anode  and  cathode  channel 


=  0,  Hs  (solid  contacts)  At :  ad,  eh 
3T.  ar  aT  M 

— -  =/lair(T’air-Tic)U 


/lairKTair  -  Ts)wfc  +  (Pair  -  AX^fc  +  Hfc)l  -  tOfclVoj Cp,o2 7*  =  Cp, 
hf[(Tt  -  Ts)wIc  +  (Tf  -  Tic)(wtc  +  Htc)]  +  wtc  ^  N,Cp,i Tt  =  Cp,{- 
z  =  0 

Pair  =  Tf  =  1073  K 
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ative  to  the  anode.  The  gases  flowing  on  each  side  of  the  PEN 
structure  remove  heat  from  each  surface  through  convection  as 
shown  in  the  boundary  conditions  in  Eqs.  (17),  (20)  and  (21).  The 
diffusion  of  species  also  contributes  to  the  transfer  of  heat  from 
the  PEN  structure  to  the  gas  phases.  Within  the  solid  structure 
heat  transfer  occurred  only  through  conduction  as  accounted  for 
in  Eqs.  (15)  and  (18).  The  gases  enter  the  channels  at  1073  K 
(Eq.  (22))  and  all  the  sides  of  the  solid  structure  are  insulated 
(Eqs.  (16)  and  (19)). 

The  rate  of  heat  generation  Q  was  determined  from 

Q  =  RsrAHsr  +  RwgsAHwgs  +  i  ^  -  Peel i (23) 

where  Q  stands  for  the  rate  of  heat  generation,  AH  the  enthalpy 
of  the  steam  reforming  (SR),  water-gas-shift  (WGS),  or  oxi¬ 
dation  (Ox)  reactions,  i  the  current  density  at  each  point,  n  the 
number  of  electrons  transferred  through  the  reaction,  F  the  Fara¬ 
day  constant,  V  the  voltage  of  the  cell,  and  R  is  the  rate  of  steam 
reforming  or  water-gas  shift,  per  surface  area  of  the  catalyst. 

There  is  an  entrance  effect  as  the  gas  enters  the  channels  on 
either  side  of  the  PEN  structure.  In  the  energy  balance  analy¬ 
sis,  axial  variations  of  the  convective  heat  transfer  coefficient, 
h,  therefore,  are  determined  and  taken  into  account  [33].  The 
variation  in  the  value  of  h  is  larger  on  the  cathode  side  because 
the  Reynolds  number  is  higher  for  air. 

2.3.  Electrochemistry 

Electrochemical  calculations  were  used  to  estimate  current 
densities  throughout  the  cell.  The  reversible  cell  voltage  (Eoev) 
or  the  “open  circuit  voltage”  (OCV)  is  a  function  of  temperature 
and  reactant  concentrations.  The  difference  between  the  OCV 
and  the  cell  operating  voltage  (Vceii)  equals  the  sum  of  ohmic 
losses  (the  current  density  (i)  times  the  cell  resistance  (R0hm)), 
the  activation  losses  or  overpotential  (>]aci),  and  the  concentration 
overpotentials  (rjconc)  as  follows: 

EoCV  -  Vcell  =  Rohmi  +  rjact  +  Oconc  (24) 

The  concentration  overpotentials  can  be  relatively  large 
(>0.02  V)  when  mass  transfer  is  the  rate-limiting  step  of  the 
electrochemical  reaction  within  the  cell.  For  the  present  case, 
however,  Ht  is  readily  available  throughout  the  anode  and  the 
maximum  current  density  is  sufficiently  low  that  the  concentra¬ 
tion  overpotential  was  negligible.  This  assumption  is  justified 
using  the  results  from  Aguiar  et  al.  [25]  and  also  by  perform¬ 
ing  calculations  at  areas  having  high  current  densities  and  low 
hydrogen  concentrations.  Thus,  for  the  present  model  only  acti¬ 
vation  and  ohmic  losses  were  considered. 

The  open  circuit  voltage,  Eocv,  was  evaluated  solving  the 
Nemst  equation  as  follows: 

*ocv  =  *»+^>n(^)  <25, 

where  E°  represents  the  theoretical  potential  of  the  reaction  at 
temperature  T  and  partial  pressures  /jh2,  Eo2,  Pu2o  for  hydro¬ 
gen,  oxygen  and  water,  respectively,  and  R  is  the  gas  constant. 


Table  4 

Equations  used  for  calculations  of  activation  overpotentials  [1] 


The  equations  used  for  exchange  current  densities  and 
overpotential  calculations  are  fisted  in  Table  4.  For  activa¬ 
tion  overpotential  calculations,  the  Butler- Volmer  equation  was 
applied  separately  to  the  anode  and  the  cathode.  The  values  of  the 
activation  overpotential  depend  on  the  nature  of  the  reaction,  the 
catalyst  and  the  temperature.  These  dependences  are  included 
in  the  values  of  exchange  current  densities,  i'o.  The  data  used  for 
the  overpotential  calculations  were  obtained  from  the  work  of 
Costamagna  and  co-workers  [1,34], 

3.  Results  and  discussion 

In  this  section,  the  model  results  for  the  concentration  and 
current  density  profiles  along  the  length  of  the  cell  will  be  dis¬ 
cussed  first.  The  results  for  the  methane  to  steam  ratio  will  then 
be  presented  and  the  susceptibility  of  different  areas  of  the  cell 
to  deposition  of  carbon  based  on  this  ratio  will  be  examined. 
The  temperature  distributions  along  the  length  and  through  the 
cross-section  of  the  unit  cell  will  be  given  and  the  resulting  tem¬ 
perature  gradients  shown.  Finally  the  effects  of  recycling  of  the 
anode  exhaust  gas  on  both  the  temperature  distribution  and  the 
methane  to  steam  ratio  will  be  discussed. 

3.1.  Concentration  and  current  density  profile 

Average  gas  phase  concentrations  of  reactants  and  products 
along  the  length  of  the  cell  are  shown  in  Fig.  2.  The  rapid 
depletion  of  CH4  and  H2O  and  the  corresponding  increase  in 


Fig.  2.  Bulk-mean  anode  gas  phase  composition  along  the  length  of  the  cell.  In 
this  and  all  following  figures,  the  relevant  section  of  the  unit  cell  is  indicated  in 
the  inset. 
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z  -  Dimensionless  Distance 


Fig.  3.  Current  density  in  the  mid-plane  of  the  PEN  structure  along  the  length 
of  the  cell. 


z  -  Dimensionless  Distance 


Fig.  4.  Ratio  of  methane  to  steam  in  the  gas  phase  along  the  length  of  the  cell. 


H2  concentration  indicates  that  steam  reforming  occurs  mainly 
at  the  entrance  to  the  cell  (z<  0.2).  The  produced  H2  is  then 
consumed  by  the  electrochemical  reaction  such  that  the  H2  con¬ 
centration  reaches  a  maximum  at  approximately  20%  along  the 
length  of  the  cell  (z  =  0.2).  Correspondingly,  there  is  a  mini¬ 
mum  in  the  H2O  concentration  because  H2O  is  consumed  by 
the  reforming  reaction  and  produced  by  the  H2  oxidation  reac¬ 
tion.  The  distributions  of  CO  and  CO2  along  the  length  of  the 
cell  are  in  agreement  with  equilibrium  calculations  for  z  >  0.2.  At 
the  inlet  region  (z  <  0.1),  however,  CO  is  produced  so  quickly  by 
steam  reforming  that  the  water-gas  shift  reaction  does  not  reach 
equilibrium.  These  trends  in  concentration  are  in  agreement  with 
the  calculations  of  Aguiar  et  al.  [25]. 

The  current  densities,  calculated  in  the  mid-plane  of  the  PEN 
structure,  are  shown  in  Fig.  3.  The  current  density  increases  to 
a  maximum  at  z  «  0.25.  The  maximum  in  the  current  density 
does  not  occur  at  the  same  location  as  the  maximum  in  the  H2 
concentration  because  the  current  density  depends  on  the  tem¬ 
perature  as  well  as  the  concentration  of  H2  (see  equations  in 
Table  4).  That  is,  the  temperature  of  the  anode  changes  along 
the  length  of  the  cell  as  the  relative  rates  of  the  exothermic  and 
endothermic  reactions  change.  The  temperature  profile  will  be 
discussed  in  Section  3.3. 

3.2.  Methane  to  steam  ratio 

In  terms  of  carbon  deposition,  the  main  parameters  of  interest 
are  the  concentration  of  methane  and  the  steam  to  hydrocarbon 
ratio.  The  steam  to  hydrocarbon  ratio  increases  along  the  length 
of  the  cell  as  the  steam  and  methane  react  to  form  hydrogen, 
which  then  electrochemically  reacts  to  form  steam  (Fig.  2).  In 
Fig.  4,  the  data  are  presented  in  terms  of  the  methane  to  steam 
ratio,  which  varies  between  0  and  0.5  (the  steam  to  methane  ratio 
varies  between  infinity  and  2).  At  the  inlet  of  the  cell,  the  ratio  of 
methane  to  steam  is  0.5,  and  this  ratio  decreases  throughout  the 
length  of  the  cell  as  methane  is  consumed  and  steam  is  produced. 
Thus,  the  inlet  of  the  cell  is  the  most  susceptible  area  for  carbon 
formation. 


The  cross-sectional  distribution  of  the  methane  to  steam  ratio 
within  the  anode  at  the  inlet  (z  =  0.005)  is  illustrated  in  the  con¬ 
tour  plot  in  Fig.  5.  The  letters  (a,  d,  e,  and  h)  correspond  to  the 
same  letters  used  in  Fig.  lb.  The  rate  of  change  in  the  methane  to 
steam  ratio  with  distance  into  the  anode  indicates  that  the  major¬ 
ity  of  the  steam  reforming  reaction  occurs  in  the  top  portion  of 
the  anode,  closest  to  the  fuel  channel  and  furthest  from  the  triple¬ 
phase  boundary,  at  which  the  fuel,  oxygen  anions  and  catalyst 
meet.  The  methane  to  steam  ratio  continuously  decreases  from 
the  gas/anode  interface  (line  ad)  through  the  anode  structure  to 
the  anode/electrolyte  interface  (close  to  line  eh — the  electrolyte 
and  cathode  are  very  thin  relative  to  the  anode).  According  to  the 
methane  to  steam  ratio,  the  area  that  is  most  susceptible  to  car¬ 
bon  formation  is  the  catalyst  at  the  gas/anode  interface  at  which 
the  steam  to  carbon  ratio  is  lowest,  and  to  which  the  oxygen 
anions  will  not  penetrate. 


Fig.  5.  Contour  plot  of  methane  to  steam  ratio  in  the  XY  plane  of  the  anode 
at  z=  0.005.  Note  that  the  values  for  x  and  y  are  the  actual  distances.  The 
concentration  of  methane  is  essentially  zero  under  the  interconnect. 
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z  -  Dimensionless  Distance 

Fig.  6.  Bulk-mean  temperature  in  the  gas  phases  along  the  length  of  the  cell 
and  temperature  along  the  mid-plane  of  the  PEN  structure.  Because  the  PEN 
structure  is  thin,  temperature  variation  in  the  y-direction  is  negligible. 

3.3.  Temperature  distribution 

In  this  section,  temperature  distributions  in  three  dimensions 
for  each  of  the  gas  and  solid  phases  are  discussed.  Calculated 
bulk-mean  gas  phase  temperature  distributions  are  shown  in 
Fig.  6.  Also  included  is  the  temperature  profile  in  the  flow  direc¬ 
tion  at  the  mid-plane  of  the  solid  phase  (PEN  structure).  The 
temperature  of  the  solid  phase  is  lowest  at  the  inlet  where  the  rate 
of  the  endothermic  steam  reforming  reaction  is  the  highest.  The 
solid  phase  is  approximately  50  K  cooler  than  the  temperature 
of  the  inlet  gases.  As  the  exothermic  electrochemical  oxida¬ 
tion  reaction  rate  increases,  the  temperature  increases  along  the 
length  of  the  cell.  The  temperature  of  the  gas-phase  on  the  anode 
side  is  close  to  the  temperature  of  the  solid  phase.  There  is  a  lag 
in  the  temperature  at  the  gas  inlet  while  the  gas  cools  to  the 
temperature  of  the  solid.  On  the  cathode  side,  the  gas  flow  rate 
is  substantially  higher  than  that  on  the  anode  side.  Changes  in 
the  air  temperature  are  therefore  less  dramatic  than  that  for  the 
anode  gas  phase.  Temperature  gradients  in  the  PEN  structure, 
along  the  length  of  the  cell,  are  relatively  large  and  attain  a 
maximum  value  of  16.4  K  cm- 1  at  z  ^  0.25. 

A  balance  on  the  total  heat  required,  taking  into  account  both 
the  endothermic  reforming  reactions  and  the  exothermic  elec¬ 
trochemical  and  WGS  reactions,  indicates  that  even  if  the  total 
amount  of  cooling  could  be  distributed  evenly  throughout  the 
cell,  the  outlet  temperature  would  be  approximately  1120K. 
That  is,  a  net  excess  of  heat  is  generated.  The  maximum  temper¬ 
ature  gradient,  however,  would  decrease  to  ~7-9  K  cm- 1  as  the 
minimum  temperature  would  be  1073  K,  the  inlet  temperature 
of  the  gases,  compared  to  1023  K  in  the  current  configuration. 

Both  the  longitudinal  (z-axis)  distribution  of  temperature  and 
the  cross-sectional  profiles  (x-y  plane)  of  temperature  determine 
thermal  stresses  in  the  materials.  Calculations  show  that  the  max¬ 
imum  temperature  difference  across  the  PEN  structure  (in  the 
y-direction)  is  within  1  K.  Therefore,  gradients  in  the  y-direction 
are  small,  as  would  be  expected  for  a  solid  with  high  thermal 
conductivity  and  with  a  relatively  thin  PEN  structure.  Heat  is 
also  well  dissipated  through  the  interconnects. 


Fig.  7.  Temperatures  in  the  PEN  structure  vs.  position  in  the  direction  lateral  to 
that  of  flow.  Z  represents  the  dimensionless  distance  from  the  gas  inlet,  along 
the  flow  direction.  Temperature  variations  in  the  y-direction  are  negligible. 


Along  the  length  of  the  cell,  however,  the  distribution  of 
temperature  within  the  PEN  structure  in  the  x-di  recti  on  varies 
considerably.  To  illustrate  the  gradients  in  the  x-direction,  the 
temperature  distributions  of  four  planes  from  z  =  0.005  (inlet 
of  cell)  to  z  =  0.5  (midway  along  length  of  cell)  are  plotted  in 
Fig.  7.  At  the  inlet,  the  temperature  difference  between  the  PEN 
structure  under  the  interconnect  and  the  PEN  structure  at  the 
middle  of  the  cross-section  is  10  K  with  a  maximum  temper¬ 
ature  gradient  of  35  Kcm-1.  Mid- way  along  the  length  of  the 
cell,  however,  this  temperature  difference  is  less  than  1 K  and 
the  temperature  gradient  is  negligible.  The  temperature  gradi¬ 
ents  in  the  x-direction  are  a  result  of  the  fact  that  only  part  of 
the  PEN  structure  is  in  contact  with  the  interconnect,  which 
effectively  dissipates  heat.  Note  that  the  maximum  values  of  the 
temperature  gradients  in  the  x-direction  are  significantly  larger 
than  those  in  the  z-direction  along  the  pathway  of  the  gases  (i.e., 
35  Kcm-1  versus  16  Kcm-1).  This  fact  illustrates  the  signifi¬ 
cance  and  importance  of  3D  modeling  in  systems  that  include 
internal  reforming. 

These  calculations  were  performed  using  a  co-current  flow  of 
gases.  Running  the  model  using  a  counter-current  flow  pattern 
resulted  in  significantly  higher  temperature  gradients,  especially 
at  the  inlet  region  of  the  cell.  The  maximum  value  of  the  gradients 
in  the  counter-current  configuration  was  more  than  40%  higher 
than  in  a  co-current  flow  system.  These  results  are  in  agreement 
with  those  reported  by  Aguiar  et  al.  [25]  for  planar  systems  and 
by  Li  and  Chyu  [35]  for  tubular  systems. 

3.4.  Applicability  of  results  to  stack  analysis 

For  the  heat  balance  calculations  all  the  cell  outer  surfaces 
were  assumed  to  be  insulated.  This  assumption  is  applicable  to 
the  case  of  a  single  stand  alone  cell.  For  cells  in  a  stack,  however, 
the  heat  transfer  from  one  cell  to  another  also  should  be  taken 
into  account,  as  the  top  surface  of  each  cell  is  in  contact  with 
the  bottom  surface  of  the  next  cell.  The  calculated  temperatures 
on  the  top  and  bottom  surfaces  of  the  cell  along  the  pathway 
are  plotted  in  Fig.  8.  As  shown,  these  values  are  within  5  K  (and 
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Fig.  8.  Temperature  of  the  top  and  bottom  surfaces  of  the  interconnect  plates 
along  the  length  of  the  cell. 

actually  within  1 K  for  over  90%  of  the  cell)  throughout  the  cell, 
such  that  there  will  be  little  if  any  heat  transfer  between  cells 
in  a  stack,  and  the  distribution  of  temperature  within  the  PEN 
structure  calculated  in  this  study  is  applicable  to  either  a  single 
cell  or  a  cell  within  a  stack.  To  further  verify  this  result,  a  “two- 
adjacent-cell”  system  was  modeled.  Although  the  temperature 
distribution  in  the  interconnects  changed  slightly,  the  tempera¬ 
ture  distribution  in  the  PEN  structure  for  the  two-adjacent-cell 
system  was  identical  to  the  distribution  obtained  for  the  single 
cell  system. 

3.5.  Recycling 

The  effects  of  recycling  on  the  temperature  distribution  for 
50,  60  and  75%  recycling  of  the  anode  exhaust  gas  are  shown 
in  Fig.  9.  The  percent  of  recycling  is  defined  as  the  amount  of 
recycled  gas  divided  by  the  total  amount  of  anode  outlet  gas  mul¬ 
tiplied  by  100%.  As  mentioned  in  Section  1,  recycling  results  in 
a  dilution  of  the  methane  concentration  and  enrichment  in  the 
hydrogen  concentration  in  the  feed.  The  higher  hydrogen  con¬ 
centration  in  the  feed  results  in  a  higher  current  density  at  the 
inlet  region  with  a  corresponding  increase  in  the  temperature. 
The  higher  the  fraction  of  recycling  the  smaller  the  temperature 
difference  between  the  inlet  (z  =  0)  and  outlet  (z  =  1)  of  the  unit 
cell  (Fig.  9a).  The  maximum  temperature  gradient  decreased 
from  16  Kcm-1  for  no  recycling  to  almost  10  Kcm-1  for  a 
recycling  ratio  of  75%. 

Recycling  also  influences  the  cross-sectional  temperature 
distribution  of  the  PEN  structure  (Fig.  9b).  Note  that  the  tem¬ 
perature  difference  is  plotted  in  the  y-axis  of  Fig.  9b  for  ease 
of  comparison.  The  absolute  temperature  of  the  PEN  structure 
will  be  different  for  each  level  of  recycling.  Consistent  with 
the  result  in  the  z-direction,  increasing  the  amount  of  recy¬ 
cling  decreases  the  temperature  difference  in  the  x-direction 
within  the  PEN  structure.  The  temperature  gradients  at  the  inlet 
region  (z  =  0.005)  decreased  from  35  Kcm-1  with  no  recycling 
to  11  Kcm-1  with  75%  recycling. 

The  steam  concentration  within  the  cell  is  also  a  function  of 
the  amount  of  recycling.  With  a  higher  concentration  of  hydro¬ 
gen  in  the  feed,  the  rate  of  the  electrochemical  reaction  increases 
and  more  steam  is  produced.  The  calculated  inlet  concentra¬ 
tions  and  average  current  densities  for  various  levels  of  recycling 
are  listed  in  Table  5.  With  50  and  60%  recycling,  the  steam  to 
methane  ratio  at  the  inlet  was  the  same  as  that  with  no  recy- 


Fig.  9.  Temperature  distribution  of  the  PEN  structure  for  different  levels  of 
recycling:  (a)  along  the  length  of  the  cell  and  (b)  through  the  cross-section  at 
z  =  0.005.  The  temperature  variation  in  the  v-direction  is  negligible. 


cling.  For  75%  recycling,  the  steam  to  methane  ratio  increased. 
Although  an  increase  in  steam  concentration  is  advantageous 
in  terms  of  reducing  carbon  formation,  the  relative  amount  of 
hydrogen  decreases  as  the  amount  of  steam  increases,  with  a 
resulting  decrease  in  the  average  current  density  of  the  sys¬ 
tem.  Thus,  an  intermediate  level  of  recycling  may  be  optimal 
depending  on  the  system  parameters. 

The  susceptibility  to  carbon  formation  will  be  influenced  by 
recycling  the  anode  exhaust  gases  in  several  ways.  First,  the 
hydrogen  to  hydrocarbon  ratio  is  higher  in  the  feed.  Second, 
the  increased  hydrogen  content  will  result  in  a  higher  amount 

Table  5 

Inlet  concentration  and  average  current  density  as  a  function  of  recycling  level 

Level  of  Inlet  concentration  (mol%)  Average  current 

recycling  (%)  density  (A  cm-2) 

CH4  h2o  h2  co  co2 

0  33.0  67.0  0  0  0  0.62 

50  21.7  43.5  13.0  8.5  13.3  0.66 

60  16.5  34.4  24.6  11.0  13.5  0.71 

75  10  42  18  12.5  17.5  0.55 
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of  steam  at  the  triple-phase  boundary.  Third,  the  increase  in 
the  concentration  of  CO2  in  the  feed  will  increase  the  rate  of 
carbon  removal  through  reaction  (4).  Offsetting  these  effects, 
however,  is  the  increase  in  the  amount  of  CO  with  increasing 
levels  of  recycling.  The  WGS  reaction  is  equilibrium  limited  at 
1073  K  and,  thus,  the  ratio  of  CO  to  CO2  increases  as  the  level  of 
recycling  increases.  The  increasing  amount  of  CO  will  impede 
removal  of  carbon  from  the  surface  by  shifting  the  equilibrium 
of  reaction  (2). 

4.  Conclusions 

A  3D  model  has  been  developed  to  investigate  the  ther¬ 
mal  and  electrochemical  behavior  of  an  anode-supported  SOFC 
operating  nominally  at  1073  K  with  direct  internal  reforming  of 
methane.  The  model  was  developed  for  a  single  cell  but  is  appli¬ 
cable  to  the  analysis  of  a  stack  because  of  the  minimal  (<5  K) 
temperature  difference  between  the  top  and  bottom  of  the  inter¬ 
connect  structures.  The  results  of  the  developed  model  indicate 
that  significant  temperature  gradients  exist  both  along  the  length 
of  the  cell  (up  to  16Kcm-1)  and  within  the  cross-section  of 
the  cell  (up  to  35  Kcm-1),  specifically  at  the  inlet  to  the  cell. 
The  identification  of  the  temperature  gradients  within  the  cross- 
section  of  the  cell  is  only  possible  with  the  use  of  a  3D  model. 
The  2D  models  previously  presented  in  the  literature  are  defi¬ 
cient  in  this  respect.  The  relatively  large  temperature  gradients 
are  a  result  of  the  relative  rates  of  the  endothermic  methane  steam 
reforming  reaction  and  exothermic  electrochemical  reaction  at 
different  locations  in  the  PEN  structure.  These  large  temperature 
gradients  can  be  reduced  by  recycling  a  percentage  of  the  anode 
exhaust  gases,  which  changes  the  composition  of  the  feed.  The 
increased  concentration  of  H2  increases  the  (exothermic)  elec¬ 
trochemical  reaction,  thus  reducing  the  temperature  decrease  at 
the  inlet  of  the  cell,  while  the  increased  concentrations  of  H2, 
H2O  and  CO2  may  also  decrease  the  amount  of  carbon  deposi¬ 
tion.  There  is  an  optimal  level  of  recycling,  however,  in  terms 
of  minimizing  the  temperature  gradients  and  carbon  formation, 
and  maximizing  the  average  current  density. 

Acknowledgements 

The  authors  acknowledge  funding  support  from  the  Alberta 
Energy  Research  Institute  and  Western  Economic  Diversifica¬ 
tion  for  this  project. 

References 

[1]  P.  Costamagna,  A.  Selimovic,  M.  Del  Borghi,  G.  Agnew,  Chem.  Eng.  J. 

102  (2004)  61. 


[2]  M.  Boder,  R.  Dittmeyer,  J.  Power  Sources  155  (2006)  13-22. 

[3]  K.  Ahmed,  K.  Foger,  Catal.  Today  63  (2000)  479-487. 

[4]  P.  Heidebrecht,  K.S.  Sundmacher,  J.  Power  Sources  145  (2005)  40-49. 

[5]  R.  Fellows,  J.  Power  Sources  71  (1998)  281-287. 

[6]  H.  He,  J.M.  Hill,  Appl.  Catal.  A:  Gen.  317  (2007)  284. 

[7]  H.S.  Bengaard,  J.K.  Norskov,  J.  Sehested,  B.S.  Clausen,  L.P.  Nielsen,  A.M. 
Molenbroek,  J.R.  Rostrup-Nielsen,  J.  Catal.  209  (2002)  365-384. 

[8]  C.M.  Finnerty,  N.J.  Coe,  R.H.  Cunningham,  R.M.  Ormerod,  Catal.  Today 
46  (1998)  137-145. 

[9]  Y.B.  Lin,  Z.L.  Zhan,  J.  Liu,  S.A.  Barnett,  Solid  State  Ionics  176  (2005) 
1827-1835. 

[10]  K.M.  Walters,  A.M.  Dean,  H.Y.  Zhu,  R  J.  Kee,  J.  Power  Sources  123  (2003) 
182-189. 

[11]  S.L.  Douvartzides,  F.A.  Couteberis,  A.K.  Demin,  P.E.  Tsiakaras,  AIChE  J. 
49  (2003)  248-257. 

[12]  K.  Eguchi,  H.  Kojo,  T.  Takeguchi,  R.  Kikuchi,  K.  Sasaki,  Solid  State  Ionics 
152/153  (2002)  411. 

[13]  S.  Park,  R.  Craciun,  J.M.  Vohs,  RJ.  Gorte,  J.  Electrochem.  Soc.  146  (1999) 
3603-3605. 

[14]  J.H.  Koh,  B.S.  Kang,  H.C.  Lim,  Y.S.  Yoo,  Electrochem.  Solid  State  Lett.  4 
(2001)  A12. 

[15]  C.H.  Bartholomew,  Catal.  Rev.  Sci.  Eng.  24  (1982)  67-112. 

[16]  A.  Abudula,  M.  Ihara,  H.  Komiyama,  K.  Yamada,  Solid  State  Ionics  86-8 
(1996)  1203-1209. 

[17]  T.  Iida,  M.  Kawano,  T.  Matsui,  R.  Kikuchi,  K.  Eguchi,  J.  Electrochem.  Soc. 
154  (2007)  B234. 

[18]  J.  Kapicka,  N.I.  Jaeger,  G.  Schulz-Ekloff,  Appl.  Catal.  A:  Gen.  84  (1992) 
47. 

[19]  R.  Peters,  R.  Dahl,  U.  Kluttgen,  C.  Palm,  D.  Stolten,  J.  Power  Sources  106 
(2002) 238-244. 

[20]  V.  Hacker,  G.  Faleschini,  H.  Fuchs,  R.  Fankhauser,  G.  Simader,  M.  Ghaemi, 
B.  Spreitz,  K.  Friedrich,  J.  Power  Sources  71  (1998)  226. 

[21]  S.H.  Clarke,  A.L.  Dicks,  K.  Pointon,  T.A.  Smith,  A.  Swann,  Catal.  Today 
38(1997)411. 

[22]  E.P.  Murray,  T.  Tsai,  S.A.  Barnett,  Nature  400  (1999)  649. 

[23]  T.  Osaki,  T.  Mori,  React.  Kinet.  Catal.  Lett.  89  (2006)  333-339. 

[24]  L.  Petruzzi,  S.  Cocchi,  F.  Fineschi,  J.  Power  Sources  118  (2003)  96. 

[25]  P.  Aguiar,  C.S.  Adjiman,  N.P.  Brandon,  J.  Power  Sources  138  (2004) 
120-136. 

[26]  W.  Lehnert,  J.  Meusinger,  F.  Thom,  J.  Power  Sources  87  (2000)  57- 
63. 

[27]  M.A.  Khaleel,  Z.  Lin,  P.  Singh,  W.  Surdoval,  D.  Collin,  J.  Power  Sources 
130  (2004)  136. 

[28]  M.  Cimenti,  J.M.  Hill,  Proceedings  of  the  10th  International  Symposium 
on  SOFCs,  Nara,  Japan,  2007. 

[29]  P.  Holtappels,  L.G.J.D.  Haart,  U.  Stimming,  I.C.  Vinke,  M.  Mogensen,  J. 
Appl.  Electrochem.  29  (1999)  561. 

[30]  Institut  fur  Werkstoffe  der  Elektrotechnik  (IWE),  The  Numerical  Simula¬ 
tion  of  Fuel  Cells  of  the  PEMFC  and  SOFC  type  with  the  Finite  Difference 
Element  Method  (FDEM).  A  report  on:  www.iwe.uni-karlsruhe.de, 
2005. 

[31]  E.  Achenbach,  E.  Riensche,  J.  Power  Sources  52  (1994)  283. 

[32]  J.M.  Moe,  Chem.  Eng.  Progr.  58  (1962)  33. 

[33]  W.  Kays,  M.  Crawford,  B.  Wiegand,  Convective  Heat  and  Mass  Transfer, 
McGraw-Hill,  2005. 

[34]  P.  Costamagna,  K.  Honegger,  J.  Electrochem.  Soc.  145  (1998)  3995. 

[35]  P.W.  Li,  M.K.  Chyu,  J.  Power  Sources  124  (2003)  487. 


